Analítica Avanzada de Datos

Challenge #2

Autor/a

Afiliación

Jorge Borja Serrano

Universidad del Norte, Barranquilla

Introducción

En esta práctica ajustaremos y validaremos un modelo de RLM, y finalmente realizaremos predicción.

Contexto Analítico

The Tire Company (TTC) fabrica compuestos para llantas de F1. En la actualidad, TTC desarrolla un nuevo compuesto para la temporada 2023, con miras a superar nuevamente a su más desafiante competidor.

En el proceso de fabricación existen 3 compuestos \(x_1, x_2\) y \(x_3\) que conforman cada llanta. De acuerdo con las regulaciones, una llanta de F1 debe pesar no más de 13 kilogramos sin incluir el rin y los sensores. Además, se sabe que el costo por cada gramo de material \(j\), en dólares, es \(c_j = \sqrt{j}\), \(j=1,2,3\). TTC está convencida de que modificando la composición de la llanta es posible alcanzar un mejor rendimiento (medido en kilómetros recorridos), de tal forma que no es necesario realizar tantas paradas en pits durante una competencia regular. En la fase de experimentación, el equipo de analítica, donde usted trabaja, encontró que la formulación actual permite recorrer \(26 \pm 7\) kilómetros con un set de llantas. Se sabe que, en libras, \(x_1\in(5, 12)\), \(x_2\in(3, 20)\) y \(x_3\in(1, 6)\). Por lo tanto, el resto del peso de la llanta debe completarse con caucho sin procesar.

El día de ayer, usted fue asignado a este proyecto con miras a optimizar la composición de las llantas en TTC, de tal manera que el número de kilómetros recorridos por cada set de llantas sea considerablemente mayor que en la actualidad.

Lectura de Datos

Para leer los datos hacemos:

Las primeras 6 filas de la base de datos d son

## primeras 6 líneas
head(d)

Análisis exploratorio

Analicemos inicialmente la correlación entre las variables disponibles:

## matriz de correlación
par(mfrow = c(1,1), mar = c(.1, .1, .1, .1))
qgraph(cor(d), graph = "cor", layout = "spring", 
       sampleSize = nrow(d), 
       legend.cex = 1, alpha = 0.05)

Numéricamente, la matriz de correlación es

## matriz de correlación
cor(d)

            y         x1           x2           x3
y   1.0000000 0.29127317 -0.762485290  0.288541500
x1  0.2912732 1.00000000  0.040792570  0.049891558
x2 -0.7624853 0.04079257  1.000000000 -0.003470425
x3  0.2885415 0.04989156 -0.003470425  1.000000000

Estos resultados indican que la correlación entre \(y\) y \(x_1\) es 0.2912, entre \(y\) y \(x_2\) es -0.762, y entre \(y\) y \(x_3\) es 0.288. En cuanto a la correlaciones entre las variables independientes, podemos concluir que estas son pequeñas, lo cual sugiere que, efectivamente, \(x_1, x_2\) y \(x_3\) son independientes.

También podemos representar las correlaciones y la distribución de cada variable en un gráfico de dispersión/correlación:

## gráfico de dispersión/correlación
ggpairs(d) + theme_minimal()

El gráfico 3D entre \(x_1\), \(x_2\) y \(y\) sería:

fig <- plot_ly(d,
        x = ~x1, 
        y = ~x2, 
        z = ~y,
        text = ~rownames(d),
        color = '#BF382A')
fig <- fig %>% add_markers()
fig <- fig %>% layout(title = '\n y vs. (x1, x2)', 
                      scene = list(xaxis = list(title = 'x1'),
                                   yaxis = list(title = 'x2'),
                                   zaxis = list(title = 'y')))
fig

Ajuste del modelo

El modelo ajustado es:

## full MLR model
fit <- lm(y ~ ., data = d)
summary(fit)

Call:
lm(formula = y ~ ., data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.7744 -2.2090 -0.0784  1.8554  7.5166 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 25.24199    1.48841  16.959  < 2e-16 ***
x1           1.19280    0.15717   7.589 3.49e-12 ***
x2          -1.65137    0.08684 -19.016  < 2e-16 ***
x3           1.73492    0.26130   6.640 5.78e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.815 on 146 degrees of freedom
Multiple R-squared:  0.7584,    Adjusted R-squared:  0.7535 
F-statistic: 152.8 on 3 and 146 DF,  p-value: < 2.2e-16

La ecuación del modelo ajustado es

\[ \hat{y} = 25.24 + 1.193 x_1 - 1.651x_2 + 1.735x_3 \]

Inferencia para \(\mathbf{\beta}\)

Los intervalos de confianza del 95% para los coeficientes pueden obtenerse a través de la función confint.default() haciendo

## 95% CI para los coeficientes
confint.default(fit)

                 2.5 %    97.5 %
(Intercept) 22.3247667 28.159218
x1           0.8847651  1.500843
x2          -1.8215743 -1.481163
x3           1.2227821  2.247048

Multicolinealidad

En términos generales, el concepto de multicolinealidad es sinómino de redundancia en las variables independientes.

Uno de los supuestos fuertes del modelo de RLM es que las variables \(X_1, X_2,\ldots,X_n\) son independientes. Cuando esto no ocurre, los estimadores de \({\beta} = (\beta_1,\beta_2,\ldots,\beta_k)\) tienen propiedades distintas a las ya conocidas.

Desde el punto de vista formal, la existencia de multicolinealidad puede probarse utilizando el ill-condition number (ICN)

## ICN
kappa(fit)     

[1] 44.46287

Teniendo en cuenta que el ICN es \(>30\), aparentemente, existe multicolinealidad entre \(x_1, x_2\) y \(x_3\).

Ahora, si estamos interesados en determinar cuál de la(s) variable(s) independiente(s) con mayor grado de colinealidad, utilizamos el VIF:

## VIF
car::vif(fit)

      x1       x2       x3 
1.004185 1.001697 1.002526 

Al analizar el VIF, es claro que ningún valor es \(>5\).

Recientemente, se han implementado otras pruebas de multicolinealidad en el paquete mctest:

## otras pruebas de multicolinealidad
mctest(fit)$odiags

                         results detection
Determinant           0.99582063         0
Farrar Chi-Square     0.61635301         0
Red Indicator         0.03726144         0
sum of Lambda Invers  3.00840801         0
Theil Indicator      -1.50847225         0
Condition Number     15.61480762         0

De acuerdo con estos resultados, no es posible hablar de la existencia de multicolinealidad. En parrticular, no aparece el valor de 1 en la columna detection.

Validación de supuestos

La validación de los supuestos puede hacerse via valores \(p\) como se muestra a continuación.

Normalidad

## prueba de Normalidad
shapiro.test(rstudent(fit))$p.value

[1] 0.8899012

Como el valor \(p\) es \(>0.05\), entonces los errores del modelo siguen una distribución Normal.

Varianza constante

## prueba de varianza constante
car:::ncvTest(fit)$p

[1] 0.3356801

Como el valor \(p\) es \(>0.05\), podemos concluir que los errores tienen varianza constante.

Independencia

## prueba de independencia
car:::durbinWatsonTest(fit)$p

[1] 0.302

Este resultado indica que los errores del modelo ajustado son independientes.

Identificación de outliers

Para identificar outliers usamos los residuales estudentizados del modelo:

## outliers?
r <- rstudent(fit)
res <- which(r > 3 | r < -3)
nout <- ifelse(length(res) == 0, 0, length(res))

Este resultado indica que hay 0 outliers en los datos.

Datos influenciales

Para identificar este tipo de observaciones, utilizamos la distancia de Cook.

En R procedemos de la siguiente manera:

## gráfico de la distancia de Cook
plot(fit, which = 4, las = 1)

Este resultado indica que las observaciones 44, 57 y 142 podrían considerarse influenciales.

Ahora, podemos encapsular las pruebas anteriores en una función:

## cargar función
source('https://www.dropbox.com/scl/fi/1z0gsv6g4a9eowcoqs4xf/validar_supuestos.R?rlkey=ioqrt3r3sl3vh7wa7sc82gjqh&dl=1')

Esta función fue creada por el profesor Jorge Vélez, PhD 
<https://jorgeivanvelez.netlify.app/> para la validación 
de los supuestos básicos en RLM. 
 
En caso de presentar algún inconveniente, por favor 
repórtelo a <jvelezv@uninorte.edu.co> 
 
Última modificación:  Abril 15, 2024 
 

Finalmente, utilizamos la función validar_supuestos() para validar todos los supuestos a la vez e identificar outliers:

## validar supuestos
## fit es el objeto que contiene el modelo de RLM
validar_supuestos(fit)

         outliers        normalidad  homocedasticidad     independencia 
"no hay outliers"          "cumple"          "cumple"          "cumple" 

Conclusión: Los supuestos de Normalidad, independencia y varianza constante de los errores se cumplen. Por lo tanto, el modelo es válido para predecir. Además, parecen no existir outliers en los datos.

Mejora de la función validar_supuestos()

A continuación, cargaremos la función que hemos creado en nuestro repositorio de GitHub, la cual es una herramienta más completa para la validación y construcción de un modelo de regresión lineal múltiple. En ella seguiremos una serie de pasos para la construcción del modelo, omitiendo el primer paso de ajuste inicial.

1. Ajuste
2. Validación
   • Prueba de significancia global
   • Prueba de significancia marginal
   • \(R^2\) ajustado
3. Verificación de supuestos del error
   • Normalidad
   • Independencia
   • Varianza constante
4. Multicolinealidad
   • ICN
   • VIF
   • …
5. Outliers/Incluenciales
   • Distancia de Cook
   • DFBetas

Para utilizar esta función en nuestro entorno de trabajo, primero cargaremos el código desde el repositorio:

source("https://raw.githubusercontent.com/unfresh25/Challenge--2/master/function.R")

Esta función fue creada por Jorge Borja 
<https://jorgeborja-portfolio.vercel.app/> para la validación 
de la construcción de un modelo de regresión lineal múltiple. 
 
En caso de presentar algún inconveniente, por favor 
repórtelo a <jborja@uninorte.edu.co> 
 
Última modificación: Junio 05, 2024 
 

Luego, procederemos a evaluar el modelo (fit) utilizando la función creada, en este caso utilizamos plot = T para ver el gráfico de la distancia de cooks. Entonces,

validar_supuestos(fit, plot = T)

1. Validación

  Significancia Global Significancia Marginal            R2 ajustado 
       "Significativo"       "Significativas"        "Significativo" 

2. Verificación de residuales

      Normalidad    Independencia Homocedasticidad 
        "Cumple"         "Cumple"         "Cumple" 

3. Multicolinealidad

Multicolinealidad 
         "Cumple" 

4. Outliers/Influenciales

         Outliers     Influenciales 
"No hay outliers"     "44, 57, 142" 

Predicción

Si fuese de interés determinar el rendimiento de una llanta con compuestos

\[ \mathbf{x}_0= (10, 4, 6) \]

procedemos de la siguiente forma:

1. Creamos el vector de nuevas condiciones:

## x0
x0 <- data.frame(x1 = 10, x2 = 4, x3 = 6)
x0

2. Realizamos la estimación de \(\hat{y}|\mathbf{x}_0\)

## estimación
predict(fit, newdat = x0) 

       1 
40.97405 

Por lo tanto, \(\widehat{E[Y|\mathbf{x}_0]} = 40.97\).

3. Construimos intervalos de confianza y predicción. Para ello basta con agregar el argumento interval a predict():

El intervalo de confianza del 95% es

## confidence interval
predict(fit, newdat = x0, interval = 'confidence') 

       fit      lwr      upr
1 40.97405 38.88432 43.06377

Así las cosas, el promedio poblacional de la distancia recorrida cuando se utiliza la mezcla de componentes \(\mathbf{x}_0= (10, 4, 6)\) es

\[ \mu | \mathbf{x}_0 \in (38.88, 43.06) \] con una confianza del 95%.

Finalmente, el intervalo de predicción del 95% es

## prediction interval
predict(fit, newdat = x0, interval = 'prediction') 

       fit      lwr      upr
1 40.97405 35.03113 46.91697

Por lo tanto, la distancia recorrida por el próximo set de llantas en el que se utilice la mezcla de componentes \(\mathbf{x}_0= (10, 4, 6)\) será

\[ Y | \mathbf{x}_0 \in (35.03, 46.92) \]

Para más información sobre el cálculo de los intervalos de confianza y predicción a partir de un modelo de RLM ajustado, se recomienda consultar ?predict.lm.

Referencias

Para más detalles sobre el modelo de RLM se sugiere consultar el Capítulo 3 del texto Modelos de Regresión: Una aproximación práctica con R.